Frequency analysis method and apparatus, and spectrum spreading demodulation method and apparatus

ABSTRACT

A method for performing frequency analysis by executing discrete Fourier transform with a smaller capacity of a memory than in usual FFT processing, without reducing the frequency resolution. A-power-of-two N A =2 a  (where “a” is an integer) memories for discrete Fourier transform are prepared. A former-stage calculation step of taking out frequency components of an input signal having a-power-of-two N B =2 b  (where “b” is an integer and b&gt;a) data items, in a comb manner and of calculating N A  intermediate data items and a latter-stage calculation step of applying fast Fourier transform to the intermediate data items obtained in the former-stage calculation step, by using the N A  memories for discrete Fourier transform are executed 2 b−a  times with the frequency components taken out in the comb manner being changed, to execute discrete Fourier transform through 2 b−a  operations.

TECHNICAL FIELD

The present invention relates to methods and apparatuses for frequencyanalysis using discrete Fourier transform, and to demodulation methodsand apparatuses for spectrum spreading signals, such as GPS (globalpositioning system) satellite signals.

BACKGROUND ART

In a GPS system for measuring the position of a mobile body by usingartificial satellites (GPS satellites), a GPS receiver has a basicfunction of receiving signals from four or more GPS satellites,calculating the position of the receiver from the received signals, andinforming the user of the position.

The GPS receiver demodulates the signals sent from the GPS satellites toacquire the orbit data of the GPS satellites, and uses simultaneousequations to calculate the three-dimensional position of the receiverfrom the orbits and time information of the GPS satellites, and thedelay time of the received signals. The reason why four GPS satellitesfrom which signals are received are required is to remove the effect ofan error between the time used in the GPS receiver and the time used bythe satellites.

A commercially available GPS receiver receives a spectrum spreadingsignal radio wave called a C/A (clear and aquisition) code in an L1 bandfrom a GPS satellite (Navstar) to perform calculations for positionmeasurement.

The C/A code is a PN (pseudorandom noise) sequence code having atransmission-signal speed (chip rate) of 1.023 MHz and a code length of1023, such as a Gold code, and is a signal obtained by BPSK (binaryphase shift keying) modulating a carrier wave (hereinafter called acarrier) having a frequency of 1575.42 MHz by a signal obtained byspreading data of 50 bps. In this case, since the code length is 1023, aPN-sequence code is repeated in the C/A code with 1023 chips being usedas one period (therefore, one period is equal to 1 millisecond), asshown in FIG. 21(A).

The PN-sequence code in the C/A code differs in each GPS satellite. TheGPS receiver can detect in advance a PN-sequence code used by each GPSsatellite. In addition, the GPS receiver understands from a navigationmessage like that described later whether the receiver can receive asignal from each GPS satellite at the position of the receiver at thepoint of time. Therefore, for three-dimensional position measurement,for example, the GPS receiver receives radio waves which can be obtainedat its position at the point of time from four or more GPS satellites,applies inverse spectrum spreading to the radio waves, and performscalculations for position measurement to obtain its position.

As shown in FIG. 21(B), one bit of satellite-signal data is transferredin units of 20 periods of the PN-sequence code, that is, in units of 20milliseconds. In other words, the data transmission rate is 50 bps. ThePN-sequence code in one period, that is, 1023 chips, is inverted betweenwhen the corresponding bit is “1” and when the bit is “0.”

As shown in FIG. 21(C), one word is formed of 30 bits (600 milliseconds)in the GPS, and one sub-frame (six seconds) is formed of 10 words asshown in FIG. 21(D). In the first word of a sub-frame, as shown in FIG.21(E), a preamble having a fixed bit pattern is always inserted evenwhen data is updated, and data is transferred after the preamble.

Further, one frame (30 seconds) is formed of five sub-frames. Thenavigation message is transferred in units of one-frame data. Firstthree sub-frames in one-frame data includes information unique to asatellite, called ephemeris information. The information includes aparameter used for obtaining the orbit of the satellite, and the timewhen the satellite sent the signal.

All GPS satellites have an atomic clock, and use common timeinformation. The time when a GPS satellite sends a signal issynchronized with one second of the atomic clock. The PN-sequence codeof a GPS satellite is generated in synchronization with the atomicclock.

Orbit information in the ephemeris information is updated in units ofseveral hours. Until an update is performed, the same information isused. The orbit information in the ephemeris information can be storedin a memory of the GPS receiver so as to use the same informationprecisely for the several hours.

The navigation message included in the remaining two sub-frames in theone-frame data is information sent in common from all satellites, calledalmanac information. The almanac information is transferred in 25frames, and includes rough-position information of each GPS satelliteand information indicating which GPS satellite is available. The almanacinformation is updated by information sent from a terrestrial station.Until an update is performed, the same information is used. The almanacinformation can be stored in a memory of the GPS receiver so as to usethe same information for several months.

To receive a signal from a desired GPS satellite to obtain theabove-described data, the C/A code is phase-synchronized by using thesame PN-sequence code (hereinafter, a PN-sequence code is called just aPN code) prepared in the GPS receiver as the C/A code used by the GPSsatellite and to be received, to acquire the signal sent from the GPSsatellite, and inverse spectrum spreading is performed. When phasesynchronization with the C/A code is obtained and inverse spreading isperformed, each bit is detected, and a navigation message, includingtime information, can be obtained from the signal sent from the GPSsatellite.

The signal sent from the GPS satellite is captured by C/A-code phasesynchronization search. In the phase synchronization search, correlationbetween the PN code of the GPS receiver and the PN code of the signalreceived from the GPS satellite is detected, and when a correlationvalue obtained as the result of correlation detection is larger than avalue specified in advance, for example, it is determined that bothcodes are synchronized. When it is determined that they are notsynchronized, the phase of the PN code of the GPS receiver is controlledby some synchronization method to synchronize the PN code of the GPSreceiver with the PN code of the received signal.

Since a GPS satellite signal is obtained by BPSK modulating the carrierby a signal obtained by spreading data with a spreading code, asdescribed above, synchronization needs to be acquired for the carrierand the data in addition to the spreading code when the GPS receiverreceives the GPS satellite signal. The synchronization of the spreadingcode and that of the carrier cannot be performed independently.

The GPS receiver usually converts the carrier frequency of the receivedsignal to an intermediate frequency several megahertz from the carrierfrequency, and performs the above-described synchronization detectionprocess with a signal having the intermediate frequency. The carrier inthe intermediate-frequency signal mainly includes a frequency errorcaused by a Doppler shift corresponding to the moving speed of the GPSsatellite and an error in the frequency of a local oscillator generatedinside the GPS receiver when the received signal is converted to theintermediate-frequency signal.

Therefore, the carrier frequency in the intermediate-frequency signal isunknown due to these frequency-error factors, and the carrier frequencyneeds to be searched for. A synchronization point (synchronized phase)in one period of the spreading code depends on the positionalrelationship between the GPS receiver and the GPS satellite. Since thispositional relationship is also unknown, some synchronization method isrequired, as described above.

Conventional GPS receivers use frequency search for the carrier and asynchronization method which uses a sliding correlator, a DLL (delaylocked loop), and a costas loop.

When the above-described conventional method is used as asynchronization detection method, however, it is, in principle, notsuited to high-speed synchronization, and to compensate this, an actualreceiver needs to have multiple channels and to search for asynchronization point in parallel. When the synchronization of thespreading code and that of the carrier require time as described above,the GPS receiver has a slow response and causes inconvenience in use.

In contrast, the improvement of the capability of hardware, typical ofwhich is DSPs (digital signal processors), has implemented a method forperforming code synchronization at a high speed by using a digitalmatched filter which uses fast Fourier transform (hereinafter calledFFT) without using a sliding correlation method like that describedabove.

FFT is a version of discrete Fourier transform (DFT), implemented by ahigher-speed algorithm. Whereas discrete Fourier transform transforms asignal which changed in time to a signal in a frequency domain to allowfrequency analysis, the FFT calculation algorithm requires the samenumber (memory capacity) of memories (random access memories, RAMs) asthat of data items. When the number of data items is the 16-th power oftwo, the required capacity of a calculation RAM is about 524 KB in32-bit calculations with the real and imaginary parts of complex numbersbeing included. In the present specification, the number of memories(RAMs) means the number of memory areas each used for reading andwriting one data item, and corresponds to the capacity of a memory.

Since calculations are recursively performed while an input-data memoryis rewritten during the calculations in the FFT calculation algorithm,output data is written into the input-data memory and the input data isnot left therein. Therefore, to leave the input data, it is necessary toseparately provide the input-data memory and a calculation memory, andthe required RAM is doubled.

Even when it is not necessary to leave the input data, if the input datais a binary signal, for example, the required memory capacity for theinput data is about 65.5 KB in the above case. Since the result of FFTis not binary and the number of bits almost equal is required, a RAMhaving a capacity of about 524 KB, which is larger than the amount ofthe input data, is needed.

In frequency analysis using usual FFT calculations, the larger thenumber of data items is, the higher the frequency resolution is.However, the required capacity of a RAM is increased as described above.Since the capacity of the RAM greatly affects the chip area of an LSI(large scale integrated circuit), an increase in the capacity of the RAMcauses an increase in cost when a DSP is integrated into an LSI.Therefore, a method for maintaining the frequency resolution with asmall capacity of the RAM even if the calculation time is extended tosome extent is demanded.

In consideration of the above points, an object of the present inventionis to allow frequency analysis to be performed without reducing thefrequency resolution, by performing discrete Fourier transform with asmaller capacity of a memory than when usual FFT processing is used.

DISCLOSURE OF INVENTION

To solve the foregoing issues, a frequency analysis method according tothe present invention is

a method using a-power-of-two N_(A)=2^(a) (where “a” is an integer)memories for discrete Fourier transform, characterized in that

a former-stage calculation step of taking out frequency components of aninput signal having a-power-of-two N_(B)=2^(b) (where “b” is an integerand b>a) data items, in a comb manner and of calculating N_(A)intermediate data items; and

a latter-stage calculation step of applying fast Fourier transform tothe intermediate data items obtained in the former-stage calculationstep, by using the N_(A) memories for discrete Fourier transform

are executed a predetermined number of times with the frequencycomponents taken out in the comb manner being changed, to executediscrete Fourier transform through the predetermined number ofoperations.

An aspect of the present invention having the above-described structureuses a structure in which an FFT calculation process for the numberN_(B) of data items includes 2^(b−a) sets of FFT calculations for thenumber N_(A) of data items. In other words, when frequency components ofN_(B) input data items are taken out in a comb manner so as not tooverlap with each other, the latter-stage calculation processing for thetaken-out results can be formed of 2^(b−a) sets of exactly the same FFTcalculation structure (FFT calculations for the number N_(A) of dataitems).

Therefore, frequency components of N_(B) input data items can beanalyzed with N_(A)=2^(a) (<N_(B)) calculation memories. This means thatdiscrete Fourier transform calculations can be performed with a smallernumber of memories than in usual FFT, and frequency components can beanalyzed without reducing the frequency resolution.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram showing a frequency analysis apparatusaccording to a first embodiment of the present invention.

FIG. 2 is a view showing an FFT signal flow, used for describing afrequency analysis method according to the first embodiment of thepresent invention.

FIG. 3 is the former part of a view showing an FFT signal flow, used fordescribing the frequency analysis method according to the firstembodiment of the present invention.

FIG. 4 is the latter part of the view showing an FFT signal flow, usedfor describing the frequency analysis method according to the firstembodiment of the present invention.

FIG. 5 is a flowchart used for describing the frequency analysis methodaccording to the first embodiment of the present invention.

FIG. 6 is a block diagram showing another frequency analysis apparatusaccording to the first embodiment of the present invention.

FIG. 7 is a view showing an FFT signal flow, used for describing afrequency analysis method according to a second embodiment of thepresent invention.

FIG. 8 is the former part of a view showing an FFT signal flow, used fordescribing the frequency analysis method according to the secondembodiment of the present invention.

FIG. 9 is the latter part of the view showing an FFT signal flow, usedfor describing the frequency analysis method according to the secondembodiment of the present invention.

FIG. 10 is a view used for describing the frequency analysis methodaccording to the second embodiment of the present invention.

FIG. 11 is a block diagram showing a frequency analysis apparatusaccording to the second embodiment of the present invention.

FIG. 12 is a flowchart used for describing the frequency analysis methodaccording to the first embodiment of the present invention.

FIG. 13 is a block diagram showing another frequency analysis apparatusaccording to the first embodiment of the present invention.

FIG. 14 is a block diagram showing the structure of a GPS receiverserving as a spectrum spreading-signal demodulation apparatus accordingto an embodiment of the present invention.

FIG. 15 is a view showing an example spectrum of a correlation detectionoutput.

FIG. 16 is a view showing a general example of a method for acquiringthe synchronization between the carrier of a received signal and aspreading code.

FIG. 17 is a view showing a method for acquiring the synchronizationbetween the carrier of a received signal and a spreading code in anembodiment of the present invention.

FIG. 18 is a view showing the structure of a main section in the firstembodiment with its operation being taken into account.

FIG. 19 is a view used for describing an embodiment of the presentinvention.

FIG. 20 is a view used for describing an embodiment of the presentinvention.

FIG. 21 is a view showing the structure of a signal sent from a GPSsatellite.

BEST MODE FOR CARRYING OUT THE INVENTION

A case in which a frequency analysis method and aspectrum-spreading-signal demodulation method according to the presentinvention are applied to a GPS receiver is taken as an embodiment andwill be described below by referring to the drawings.

Frequency Analysis Method According to First Embodiment

An outline and the operation principle of a frequency analysis methodaccording to a first embodiment of the present invention will bedescribed first. In the following description, it is assumed that aninput signal having the number N_(B) of data items is Fouriertransformed for frequency analysis by using the number N_(A) ofcalculation memories (RAMs). In the present specification, the number ofmemories (RAMs) means the number of memory areas each used for readingand writing one data item.

The input signal having the number N_(B) of data items is called r(n),the result obtained by applying discrete Fourier Transform (DFT) theretois called R(k), “n” indicates a discrete time, and “k” indicates adiscrete frequency, where “n” and “k” range from zero to N_(B)−1, andN_(B) is 2^(b) (“b” is an integer). The number N_(A) of calculationmemories which can be used in a discrete Fourier transform process isset to 2^(a) (“a” is an integer).

To calculate the result R(k) of discrete Fourier transform of the inputsignal r(n) by usual fast Fourier transform (FFT), at least N_(A)=N_(B)RAMs are required to store halfway calculation data and the finalresult. In other words, input signals r(0) to r(N_(B)−1) are stored inN_(B) RAMs, the contents of the RAMs are rewritten in the process of FFTcalculations, and the final RAM values indicate R(0) to R(N_(B)−1).

In this case, since input-signal data is lost, it is necessary toseparately provide RAMs for input-signal storage and those for FFTcalculations to save the input signal. There is no special problem inthe required capacity of the RAMs when the number N_(B) of data items inthe result R(n) of discrete Fourier transform is small.

When the number N_(B) of data items increases, however, the capacity ofthe RAMs also increases and the chip area required for the correspondingLSI is also increases. This is undesirable. In addition, since the inputdata r(n) and the result R(k) of discrete Fourier transform do notnecessarily have the same signal level (number of bits) as is typicallyin a case when the input data r(n) is a binary signal, there is a casein which it is realistic that RAMs for storing the input data r(n) andRAMs for FFT calculations are separately provided.

It is not necessary to leave all of the result R(k) of discrete Fouriertransform if an object is to detect an unspecific frequency peak. Inthis case, when a method is used in which R(k) is calculated from r(n)according to the definition expression of discrete Fourier transform,and the result is sequentially checked for k=0 to N_(B)=1, sincesum-of-products calculations are performed while input data is read intoa general-purpose register, calculation RAMs are not required. In thiscase, however, since the high-speed FFT algorithm cannot be used, andthe amount of calculations becomes extraordinary large. If the number ofdata items is large, this method is unrealistic in terms of acalculation time.

In the embodiment described below, a method is used in whichdiscrete-Fourier-transform calculations are performed with thehigh-speed FFT algorithm being left and the required capacity of RAMsbeing reduced. To reduce the required capacity of RAMs to a level ofN_(A)<N_(B), that is, a<b, calculations for obtaining R(k) from r(n) areperformed in a plurality of stages as in the following procedure.

In this case, there are separately provided RAMs for storing the inputsignal r(n) and calculation RAMs for the result R(k) of discrete Fouriertransform. To divide calculations for obtaining the result R(k) ofdiscrete Fourier transform into calculations in a plurality of stages, apoint that an FFT calculation process for the number N_(B) of data itemsincludes 2^(b−a) sets of FFT calculations for the number N_(A) of dataitems is used. An example will be shown below.

FIG. 2 is a view showing an FFT signal flow in a case in which an inputsignal has eight (N_(A)=8, a=3) data items. There are shown inputsignals r(0) to r(7) at the left end and the results R(0) to R(7) ofdiscrete Fourier transform at the right end.

FIG. 3 and FIG. 4 are views showing an FFT signal flow in a case inwhich an input signal has 32 (N_(B=)32, b=5) data items. FIG. 3 showsthe first-half part of the signal flow view, and there are shown inputsignals r(0) to r(31) at the left end and intermediate calculationvalues h0(0) to h0(7), h1(0) to h1(7), h2(0) to h2(7), and h3(0) toh3(7) at the right end. FIG. 4 shows the second-half part, and there areshown the intermediate calculation values h0(0) to h0(7), h1(0) toh1(7), h2(0) to h2(7), and h3(0) to h3(7) at the left end and theresults R(0) to R(31) of discrete Fourier transform at the right end.

When FIG. 2 is compared with FIG. 4, it is understood that each of foursignal-flow portions C0, C1, C2, and C3 enclosed by thin lines in FIG. 4has the same form as the signal-flow view of FIG. 2. More specifically,the four signal-flow portions C0, C1, C2, and C3 respectively have theintermediate calculation values h0(0) to h0(7), h1(0) to h1(7), h2(0) toh2(7), and h3(0) to h3(7) as their inputs, and R0(0) to R0 (7), R1(0) toR1(7), R2(0) to R2(7), and R3(0) to R3(7) as their outputs, and the samecalculation procedure as that shown in FIG. 2 is applied in the foursignal-flow portions C0, C1, C2, and C3 although their inputs andoutputs differ from each other.

As shown in FIG. 4, C0, C1, C2, and C3, and h0(0) to h0(7), h1(0) toh1(7), h2(0) to h2(7), and h3(0) to h3(7) can be expressed by Cp andhp(0) to hp(7), where “p” (in the case shown in FIG. 4, “p”=0 in C0,“p”=1 in C1, “p”=2 in C2, and “p”=3 in C3) is the minimum value of “k”in R(k) included in the outputs R0(0) to R0(7), R1(0) to R1(7), R2(0) toR2(7), and R3(0) to R3(7) of the signal-flow portions C0, C1, C2, andC3.

A coefficient W shown in FIG. 2 is expressed by:W=exp(−j2π/8)  (expression a)and a coefficient W shown in FIG. 3 and FIG. 4 is expressed by:W=exp(−j2π/32)  (expression b)according to the number of input data items for each case. Therefore,calculations performed in FIG. 3 and those performed in the latterstages C0, C1, C2, and C3 of FIG. 4 are exactly the same. In the aboveexpressions, “j” in “j2π” indicates the imaginary unit.

As shown in FIG. 4, “k” increases from the minimum value “p” at aninterval of four in the result R(k) of discrete Fourier transformincluded in the outputs R0(0) to R1(7), R1(0) to R1(7), R2(0) to R2(7),and R3(0) to R3(7) of the signal-flow portions C0, C1, C2, and C3.

In other words, the result R(k) of discrete Fourier transform includedin the outputs R0(0) to R0(7), R1(0) to R1(7), R2(0) to R2(7), and R3(0)to R3(7) of the signal-flow portions C0, C1, C2, and C3 is arranged at aconstant interval of kd=4.

Therefore, the result R(k) of discrete Fourier transform included in theoutputs R0(0) to R0(7), R1(0) to R1(7), R2(0) to R2(7), and R3(0) toR3(7) of the signal-flow portions C0, C1, C2, and C3 includes frequencycomponents taken out in a comb manner from the frequency components ofr(n) having 32 input data items so as not to overlap with each other.

This is a characteristic derived from a point that FFT calculationsincludes multiplications of a complex sine wave and additions orsubtractions of a delayed signal in each stage. From FIG. 4, hp(0) tohp(7) are expressed by expression (1) in FIG. 19. Expression (1) of FIG.19, “i” in hp(i) is 0, 1, 2, . . . , and 7, and W=exp(−j2π/32). Theexpression (1) can be expressed collectively by expression (2) in FIG.19.

In this case, When W⁰ to W¹⁵ are stored in a ROM (read only memory) asin usual FFT, calculations related to W can be performed by using therelationships of W¹⁶=−1 and W³²=1.

According to the above description, in the first embodiment,calculations for discrete Fourier transform can be applied to input datar(n) with a memory (RAM) having a smaller capacity than in usual FFT byusing a frequency analysis apparatus having a structure shown in FIG. 1.Frequency components can be analyzed without reducing the frequencyresolution.

Specifically, the frequency analysis apparatus shown in FIG. 1 is formedof former-stage calculation means 1 for calculating an intermediatecalculation value hp(i) (in this case, h0(0) to h0(7), h1(0) to h1(7),h2(0) to h2(7), and h3(0) to h3(7)) from input data r(n) (in this case,r(0) to r(31)), latter-stage calculation means 2 provided with RAMs usedfor FFT calculations applied to the intermediate calculation valuehp(i), a calculation control section 3, and a calculation result memory4.

Under the control of the calculation control section 3, the former-stagecalculation means 1 reads input data r(i+8n) into a general-purposeregister, multiplies it by a W term, accumulates the product in anaccumulator to calculate the intermediate calculation values hp(0) tohp(7), as shown in FIG. 3.

The eight intermediate calculation value data items, which are theresults of calculations, are written into eight data calculation RAMsconstituting the latter-stage calculation means 2. The eight data itemswritten into the RAMs are converted to Rp(0) to Rp(7) by the FFTalgorithm under the control of the calculation control section 3.

In FIG. 1, R0(0) to R0(7) corresponding to p=0 are first obtained,checked, and a peak larger than a threshold, for example, is stored inthe memory 4 or output to the outside. Then, “p” is changed to one, two,and three, calculations are sequentially performed to obtain R1(0) toR1(7), R2(0) to R2(7), and R3(0) to R3(7) by the repetition ofcalculations in the same calculation RAMs for eight data items, andfrequency components thereof are sequentially checked. With this, 32frequency components R(0) to R(31) can be finally analyzed.

The foregoing description corresponds to a case in which N_(A)=8, a=3,N_(B)=32, and b=5. Expression (2) in the foregoing calculationprocessing can be generalized to obtain expression (3) shown in FIG. 19.

Since hp(i) in expression (2) and expression (3) indicates the signalobtained by applying inverse discrete Fourier transform to a signalhaving discrete frequency components (in the above case, at an intervalof four), such as Rp(0) to Rp(7) shown in FIG. 3 and FIG. 4, hp(i)corresponds to a time waveform obtained when passing through a combfilter.

In the same way as in the above-described case, when the results Rp(0)to Rp(2^(a)−1) of calculations performed 2^(b−a) times with “p” beingset to zero to 2^(b−a)−1 are sequentially checked in the discreteFourier transform calculations generalized in this way, the same numberof calculation RAMs as the number of data items N_(A)=2^(a) can be usedto analyze N_(B)=2^(b) frequency components R(0) to R(2^(b)−1), whereN_(B) is larger than N_(A).

FIG. 5 shows an example flowchart of the frequency analysis methoddescribed above.

The values of N_(A), N_(B), “a,” and “b” are first specified as initialsettings (step S101). Then, input data r(n) is input. In the case shownin FIG. 1, the input data r(n) is input to the former-stage calculationmeans 1 (step S102).

Next, “p” is set to 0 to specify the first signal-flow portion C0 amongsignal-flow portions Cp (p=0 to 2^(b−a)−1) where frequency componentsare taken out in a comb manner (step S103). Then, a general-purposeregister is used to obtain an intermediate calculation value hp(i)according to the above-described expression (3) (step S104).

Next, the obtained intermediate calculation value hp(i) is written intothe same number of calculation RAMs as the number of data itemsN_(A)=2^(a), and FFT calculations are performed to obtain Rp(k) (stepS105). The obtained Rp(k) is frequency-analyzed (step S106), and theresult of analysis is written into a memory or output to the outside(step S107).

Next, it is determined whether the variable “p” is smaller than itsmaximum value 2^(b−a)1 (step S108). When the variable “p” is smallerthan 2^(b−a)−1, the variable “p” is incremented by one (step S109), theprocedure returns to step S104, and the processes of step S104 andsubsequent steps are repeated.

With the above-described method, a smaller number of calculation RAMsthan the number of input data items are used for discrete Fouriertransform calculations to analyze the same number of frequencycomponents as the number of the input data items. In a case whereN_(A)=4,096, a=12, N_(B)=65,536, and b=16, the required number ofcalculation RAMs is one sixteenth that for usual FFT calculations.

When the number of times complex-number multiplications are performed iscompared between usual FFT and the above-described method according tothe present invention, calculations are performed 4.9×10⁵ times in theformer case whereas calculations are performed 1.4×10⁶ times in themethod according to the present invention, which is about three timeslarger. The method needs, however, a much smaller amount of calculationsthan that, 4.3×10⁹, required for calculations performed according to thedefinition expression of discrete Fourier transform, and thelatter-stage FFT calculation greatly makes the calculations faster.

In FIG. 1, the former calculation means 1 performs FFT calculations foreach of the signal-flow portions C0, C1, C2, and C3 to obtain eightintermediate calculation values, writes the obtain intermediatecalculation values in the RAMs of the latter-stage calculation means 2to execute FFT calculations. The structure of the circuits can bechanged as shown in FIG. 6.

In a case shown in FIG. 6, a former-stage calculation means 5 obtainsall intermediate calculation values sent to the signal-flow portions C0,C1, C2, and C3. The intermediate calculation values are sequentiallywritten into RAMs of a latter-stage calculation means 2 in units ofeight values through an input switching section 6 so as to be suited toFFT processing performed in the RAMs of the latter-stage calculationmeans 2.

Second Embodiment of Frequency Analysis Method

An outline and the operation principle of a frequency analysis methodaccording to a second embodiment of the present invention will bedescribed first. Also in this case, it is assumed that an input signalhaving the number N_(B) of data items is Fourier transformed forfrequency analysis by using the number N_(A) of calculation memories(RAMs).

As described before, the frequency analysis method according to thefirst embodiment requires an increased number of complex-numbermultiplications, about slightly less than three times that required inusual FFT. A frequency analysis method according to the secondembodiment improves this point.

Specifically, the second embodiment provides a method for performing ata higher speed comb-filter calculations and calculations using FFT bothdescribed in the first embodiment. As described later, an input signalis limited to a real signal (signal having no imaginary part, such as asignal received from a GPS satellite) in order to reduce the number oftimes the comb-filter calculations and FFT calculations are performed.

Also in the second embodiment, to divide calculations for obtaining theresult R(k) of discrete Fourier transform into calculations in aplurality of stages, a point that FFT calculations for the number N_(B)of data items include 2^(b−a) sets of FFT calculations for the numberN_(A) of data items is used. An example will be shown below.

FFT calculations for N_(A) data items according to the second embodimentare performed as shown in FIG. 7. Specifically, FIG. 7 is a view showingan FFT signal flow in a case in which an input signal has four (N_(A)=4,a=2) data items. There are shown input signals r(0) to r(3) at the leftend and the results R(0) to R(3) of discrete Fourier transform at theright end.

FIG. 8 and FIG. 9 are views showing an FFT signal flow in a case inwhich an input signal has 32 (N_(B=)32, b=5) data items, which wasdescribed before and shown in FIG. 3 and FIG. 4. In the case shown inFIG. 8 and FIG. 9 in the second embodiment, the first-half part of thesignal-flow view shown in FIG. 8 and the second-half part shown in FIG.9 differ from the case shown in FIG. 3 and FIG. 4 in the firstembodiment.

FIG. 8 shows the first-half part of the signal flow view, and there areshown input signals r(0) to r(31) at the left end and intermediatecalculation values t0(0) to t0(3), t1(0) to t1(3), t2(0) to t2(3), t3(0)to t3(3), t4(0) to t4(3), t5(0) to t5(3), t6(0) to t6(3), and t7(0) tot7(3) at the right end.

FIG. 9 shows the second-half part of the signal-flow view, and there areshown the intermediate calculation values t0(0) to t0(3), t1(0) tot1(3), t2(0) to t2(3), t3(0) to t3(3), t4(0) to t4(3), t5(0) to t5(3),t6(0) to t6(3), and t7(0) to t7(3) at the left end and the results R(0)to R(31) of discrete Fourier transform at the right end.

When FIG. 7 is compared with FIG. 9, it is understood that each of eightsignal-flow portions D0, D1, D2, D3, D4, D5, D6, and D7 enclosed by thinlines in FIG. 9 has the same form as the signal-flow view of FIG. 7.More specifically, the eight signal-flow portions D0, D1, D2, D3, D4,D5, D6, and D7 respectively have t0(0) to t0(3), t1(0) to t1(3), t2(0)to t2(3), t3(0) to t3(3), t4(0) to t4(3), t5(0) to t5(3), t6(0) tot6(3), and t7(0) to t7(3) as their inputs, and R0(0) to R0(3), R1(0) toR1(3), R2(0) to R2(3), R3(0) to R3(3), R4(0) to R4(3), R5(0) to R5(3),R6(0) to R6(3), and R7(0) to R7(3) as their outputs, the samecalculation procedure is applied in the eight signal-flow portions D0,D1, D2, D3, D4, D5, D6, and D7 although their inputs and outputs differfrom each other.

As shown in FIG. 9, D0, D1, D2, D3, D4, D5, D6, and D7, and t0(0) tot0(3), t1(0) to t1(3), t2(0) to t2(3), t3(0) to t3(3), t4(0) to t4(3),t5(0) to t5(3), t6(0) to t6(3 t7(0) to t7(3) can be expressed by Dq andtq(0) to tq(3), where “q” (in the case shown in FIG. 9, “q”=0 in D0,“q”=1 in D1, “q”=2 in D2, “q”=3 in D3, “q”4 in D4, “q”=5 in D5, “q”=6 inD6, and “q”=7 in D7) is the minimum value of “k” in R(k) included in theoutputs R0(0) to R0(3), R1(0) to R1(3), R2(0) to R2(3), R3(0) to R3(3),R4(0) to R4(3), R5(0) to R5(3), R6(0) to R6(3), and R7(0) to R7(3) ofthe signal-flow portions D0, D1, D2, D3, D4, D5, D6, and D7.

A coefficient W shown in FIG. 7 is expressed by:W=exp(−j2π/4)  (expression a)and a coefficient W shown in FIG. 8 and FIG. 9 is expressed by:W=exp(−j2π/32)   (expression b)according to the number of input data items for each case. Therefore,calculations performed in FIG. 7 and those performed in the latterstages D0, D1, D2, D3, D4, D5, D6, and D7 of FIG. 9 are exactly thesame. In the above expressions, “j” in “j2π” indicates the imaginaryunit.

As shown in FIG. 9, “k” increases from the minimum value “q” at aninterval of eight in the result R(k) of discrete Fourier transformincluded in the outputs R0(0) to R0(3), R1(0) to R1(3), R2(0) to R2(3),R3(0) to R3(3), R4(0) to R4(3), R5(0) to R5(3), R6(0) to R6(3), andR7(0) to R7(3) of the signal-flow portions D0, D1, D2, D3, D4, D5, D6,and D7.

In other words, the result R(k) of discrete Fourier transform includedin the outputs R0(0) to R0(3), R1(0) to R1(3), R2(0) to R2(3), R3(0) toR3(3), R4(0) to R4(3), R5(0) to R5(3), R6(0) to R6(3), and R7(0) toR7(3) of the signal-flow portions D0, D1, D2, D3, D4, D5, D6, and D7 isarranged at a constant interval of k=8.

Therefore, the result R(k) of discrete Fourier transform included in theoutputs R0(0) to R0(3), R1(0) to R1(3), R2(0) to R2(3), R3(0) to R3(3),R4(0) to R4(3), R5(0) to R5(3), R6(0) to R6(3), and R7(0) to R7(3) ofthe signal-flow portions D0, D1, D2, D3, D4, D5, D6, and D7 includesfrequency components taken out in a comb manner from the frequencycomponents of r(n) having 32 input data items so as not to overlap witheach other.

For example, FIG. 10 shows the relationship between the signal-flowportion D1 and the signal-flow portion D7. As clearly shown in FIG. 10,the result R(k) of discrete Fourier transform included in the resultsR1(0) to R1(3) of discrete Fourier, included in the signal-flow portionD1 and the results R7(0) to R7(3) of discrete Fourier, included in thesignal-flow portion D7 differs in that “k” in R(k) is apart from eachother by seven. The results do not overlap with each other, and includefrequency components taken out in a comb manner.

Therefore, also in the second embodiment, since the former-stagecalculation means obtains the intermediate value tq(m) (where m is aninteger equal to or larger than zero, and in this case, M=0, 1, 2, or 3)from the input data r(n), and the latter-stage calculation means havingRAMs used for applying FFT calculations to the intermediate calculationvalue tq(m) sequentially calculates D0 to D7, discrete Fourier transformcalculations can be applied to the input data r(n) with a smaller numberof memories (RAMs) than in usual FFT, and frequency components can beanalyzed without reducing the frequency resolution.

However, the same problem as in the above-described first embodimentoccurs in this condition. In the second embodiment, in a case in whichan input signal is limited to a real signal having no imaginary part,since the case has a characteristic described later, the latter-stagecalculation means performs FFT calculations only for the signal-flowportions D0 to D4 among the signal-flow portions D0 to D7 to obtain theresults of FFT calculations for all the signal-flow portions D0 to D7.The calculation processing speed is made faster by this technique. Amethod for making the calculation processing speed faster in the secondembodiment will be described below.

When a signal having no imaginary part is converted to a signal in thefrequency domain, it is generally known that the real part of the resultR(k) of conversion is an even function and the imaginary part is an oddfunction. In this case, the following expressions are satisfied.Rr(k)=Rr(N _(B) −k)  (expression c)Ri(k)=Ri(N _(B) −k)  (expression d)where, Rr(k) and Ri(k) indicate the real part and the imaginary part ofR(k), respectively, and N_(B) indicates the number of data items.

The signal-flow portions D1 and D7 shown in FIG. 9 are:

$\begin{matrix}{{D1} = \left\{ {{{R1}\mspace{11mu}(0)},{{R1}\mspace{11mu}(1)},{{R1}\mspace{11mu}(2)},{{R1}\mspace{11mu}(3)}} \right\}} \\{= \left\{ {{R(1)},{R(9)},{R(17)},{R(25)}} \right\}} \\{{D7} = \left\{ {{{R7}\mspace{11mu}(0)},{{R7}\mspace{11mu}(1)},{{R7}\mspace{11mu}(2)},{{R7}\mspace{11mu}(3)}} \right\}} \\{= \left\{ {{R\mspace{11mu}(7)},{R\mspace{11mu}(15)},{R\mspace{11mu}(23)},{R\mspace{11mu}(31)}} \right\}}\end{matrix}$

From the expression c, the following expression is satisfied.

$\begin{matrix}{{{R7}\mspace{11mu}(0)} = {{{R7r}\mspace{11mu}(0)} + {{jR7i}\mspace{11mu}(0)}}} \\{= {{{Rr}\mspace{11mu}(7)} + {{jRi}\mspace{11mu}(7)}}} \\{= {{{Rr}\mspace{11mu}(25)} - {{jRi}\mspace{11mu}(25)}}} \\{= {{{R1r}\mspace{11mu}(3)} - {{R1i}\mspace{11mu}(3)}}} \\{= {{R1}^{\prime}\mspace{11mu}(3)}}\end{matrix}$where, “j” is the imaginary unit, Rxr(k) and Rxi(k) indicate the realpart and the imaginary part of Rx(k), respectively, and Rx′(k) indicatesRx(k) in which the sign of the imaginary part is inverted.

In the same way, the following expressions are also satisfied.R7(1)=R1′(2)R7(2)=R1′(1)R7(3)=R1′(0)Then, the following expression is obtained.R7(m)=R40(N _(c)−1−m) (0≦m<N _(c))where, N_(c) indicates the number of data items in a signal-flowportion.

In the same way, the following relationships are satisfied between thesignal-flow portions D2 and D6, and the signal-flow portions D3 and D5.R6(m)=R2′(N _(c)−1−m)R5(m)=R3′(N _(c)−1−m)

Such a relationship is not satisfied between the signal-flow portions D0and D4. In a more generalized manner, such a relationship is notsatisfied between D0 and D[n/2] in n signal-flow portions.

As described above, when the signs of the imaginary parts of theelements of the FFT result obtained in each of some portions of thesignal-flow portions D0 to D7 are inverted and the arrangement order ofthe elements are inverted, the elements become equal to the elements ofthe FFT result obtained in another signal-flow portion. Morespecifically, when the results of FFT calculations are obtained in thesignal-flow portions D1, D2, and D3, the results of FFT calculations forthe signal-flow portions D7, D6, and D5 can be obtained by inverting thesigns of the imaginary parts of the elements of the FFT results obtainedin the signal-flow portions D1, D2, and D3 and by inverting thearrangement order of the elements.

Therefore, when the whole of the input signal is FFT processed andconverted to signals in the frequency domain, if FFT calculations areperformed only in the signal-flow portions D0 to D4 among the eightsignal-flow portions D0 to D7, the results of the remaining signal-flowportions D5 to D7 can be easily calculated from the results of FFTcalculations in the signal-flow portions D1 to D3. Consequently, FFTcalculations for the signal-flow portions D5 to D7, and comb filtercalculations used for obtaining inputs to the signal-flow portions D5 toD7 can be omitted.

With this, according to the second embodiment, FFT can be performed ahigher speed than in the first embodiment with the amount of use memoryand the frequency resolution being about equal to those in the firstembodiment. When the number of complex-number multiplications iscompared among usual FFT, the first embodiment, and the secondembodiment, in a case in which N_(A) is set to 4,096 and N_(B) is set to65,536, for example, whereas the usual FFT requires 4.9×10⁵multiplications, and the first embodiment needs 1.4×10⁶ multiplications,which is about slightly less than three times those in the usual FFT,the second embodiment just requires 7.9×10⁵ multiplications, which isabout 1.6 times those required in the usual FFT.

From the above description, with the use of a frequency analysisapparatus having a structure shown in FIG. 11, discrete Fouriertransform calculations can be applied to input data r(n) at a high speedwith a small number of memories (RAMs) as in usual FFT, and frequencycomponents can be analyzed at a high speed without reducing thefrequency resolution.

Specifically, the frequency analysis apparatus shown in FIG. 11 isformed of former-stage calculation means 11 for calculating anintermediate calculation value tq(m) (in the case shown in FIG. 8 andFIG. 9, 0≦q≦4, 0≦m≦3) from input data r(n) (in the case shown in FIG. 8and FIG. 9, r(0) to r(31)), latter-stage calculation means 12 providedwith RAMs used for FFT calculations applied to the intermediatecalculation value tq(m) (in the case shown in FIG. 8 and FIG. 9, 0≦q≦4,0≦m≦3), for inverting the arrangement orders of the elements of obtainedFFT results and inverting the signs of the imaginary parts to obtain theresults of other FFT calculations, a calculation control section 13, anda calculation result memory 14.

Under the control of the calculation control section 13, theformer-stage calculation means 11 reads input data r(n) into ageneral-purpose register, multiplies it by a W term, accumulates theproduct in an accumulator to calculate the intermediate calculationvalues tq(m), as shown in FIG. 8. The intermediate calculation valuetq(m) is not calculated for all in 0≦q≦(Q−1) (where Q indicates thenumber of signal-flow portions having N_(A) data items, and in the casedescribed above, Q=8), but in 0≦q≦Q/2. More specifically, in the caseshown in FIG. 8 and FIG. 9, five sets of intermediate calculation valuest0(0) to t0(3), t1(0) to t1(3), t2(0) to t2(3), t3(0) to t3(3), andt4(0) to t4(3) input to the signal-flow portions D0 to D4.

The (Q/2+1) sets of intermediate calculation value data, which is theresult of calculations, are written into m data calculation RAMsconstituting the latter-stage calculation means 12. The m data writteninto the RAMs are converted to Rq(m) by the FFT algorithm under thecontrol of the calculation control section 13.

In the results of calculations, R1(m) to R[Q/x−1] (m) among Rq(m), thearrangement order of the elements thereof are reversed, and the signs ofthe imaginary parts are also inverted to calculate the results of FFTcalculations, R[Q/2+1] (m) to R[Q−1] (m) in the latter-stage calculationmeans 12.

In the case described above by referring to FIG. 8 and FIG. 9, theresults of FFT calculations, R0(0) to R0(3), corresponding to q=0 arefirst obtained in FIG. 11, checked, and a peak larger than a threshold,for example, is stored in the memory 14 or output to the outside. Then,the results of FFT calculations, R1(0) to R1(3), corresponding to q=1are first obtained, checked, and a peak larger than a threshold, forexample, is stored in the memory 14 or output to the outside.

Next, the arrangement order of the elements of the results of FFTcalculations, R1(0) to R1(3), corresponding to q=1 is inverted, and thesigns of the imaginary parts of the elements are inverted to obtain theresults of FFT calculations, R7(0) to R7(3), corresponding to q=7. Theyare checked, and a peak larger than a threshold, for example, is storedin the memory 14 or output to the outside.

Next when q=2 and q=3, calculations are performed in the same way as inthe case of q=0 to obtain R2(0) to R2(3) and R3(0) to R3(3) by therepetition of calculations in the same calculation RAMs for four dataitems, and frequency components therefor are sequentially checked. Theresults are stored in the memory 14 or output to the outside.

Then, the arrangement order of the elements of each of the results ofFFT calculations, R2(0) to R2(3) and R3(0) to R3(3) is inverted, and thesigns of the imaginary parts of the elements are inverted to obtain theresults of FFT calculations, R6(0) to R6(3) and R5(0) to R5(3),corresponding to q=6 and q=5. They are checked, and a peak larger than athreshold, for example, is stored in the memory 14 or output to theoutside.

At the end, the results of FFT calculations, R4(0) to R4(3),corresponding to q=4 are obtained, checked, and a peak larger than athreshold, for example, is stored in the memory 14 or output to theoutside. With this, frequency analysis has been applied to 32 frequencycomponents R(0) to R(32).

FIG. 12 shows an example flowchart of the frequency analysis methoddescribed above. When the processing shown in the flowchart mainlycorresponds to the control processing of the calculation control section13 in the case shown in FIG. 11.

The values of N_(A), N_(B), “a,” and “b” are first specified as initialsettings (step S201). Then, input data r(n) is input. In the case shownin FIG. 11, the input data r(n) is input to the former-stage calculationmeans 11 (step S202).

Next, “q” is set to 0 to specify the first signal-flow portion D0 amongsignal-flow portions Dq (q=0 to 2^(b−a)−1) where frequency componentsare taken out in a comb manner (step S203). Then, a general-purposeregister is sued to obtain the intermediate calculation value tq(m)(step S204).

Next, the obtained intermediate calculation value tq(m) is written intothe same number of calculation RAMs as the number of data itemsN_(A)=2^(a), and FFT calculations are performed to obtain Rq(m) (stepS205). The obtained Rp(m) is frequency-analyzed (step S206).

Next, it is determined whether the variable “q” falls in a range of0<q<2^(b−a−1) (step S207). When the variable “q” does not fall in thisrange, only the result of analysis obtained in step S206 is output (stepS210). In the case described above by referring to FIG. 8 and FIG. 9, itis determined in step S207 whether one of the signal-flow portions D1 toD3, where the result of FFT for another signal-flow portion is obtainedby inverting the arrangement order of the elements of the result of FFTand inverting the imaginary parts, is handled.

Next, it is determined whether the variable “q” satisfies q≦2^(b−a−1)−1(step S211). When the variable “q” does not satisfy the expression, theprocessing routine is terminated. When the variable “q” satisfies theexpression, the variable “q” is incremented by one (step S212), theprocedure returns to step S204, and the processes of step S204 andsubsequent steps are repeated.

When it is determined in step S207 that the variable “q” falls in arange of 0<q<2^(b−a−1), the arrangement order of the elements of theresult of FFT obtained at the variable “q” is inverted and the imaginaryparts are also inverted to obtain the result of FFT obtained when thevariable “q” is (2^(b−a)−q) (step S208), and the result isfrequency-analyzed (step S209). Then, in step S210, the result ofanalysis obtained in step S206 and the result of analysis obtained instep S209 are output. The procedure proceeds to step S211.

In FIG. 11, in the case shown in FIG. 8 and FIG. 9, the formercalculation means 11 performs FFT calculations for each of thesignal-flow portions D0, D1, D2, and D3 to obtain four intermediatecalculation values, writes the obtain intermediate calculation values inthe RAMs of the latter-stage calculation means 12 to execute FFTcalculations. The structure of the circuits can be changed as shown inFIG. 13.

In a case shown in FIG. 13, in the case shown in FIG. 8 and FIG. 9, aformer-stage calculation means 15 obtains all intermediate calculationvalues t0(0) to t0(3), t1(0) to t1(3), t2(0) to t2(3), t3(0) to t3(3),and t4(0) to t4(3) sent to the signal-flow portions D0, D1, D2, D3, andD4 to which signal inputs are required. The intermediate calculationvalues are sequentially written into RAMs of a latter-stage calculationmeans 12 in units of four values through an input switching section 16so as to be suited to FFT processing performed in the RAMs of thelatter-stage calculation means 12.

Example Structure of a GPS Receiver

The frequency analysis method in the present embodiment described aboveis used for demodulation of spectrum spreading signals in a GPS receiverto make the structure of a DSP 100 simple and to allow high-speedprocessing to be expected. An example structure of a GPS receiverserving as a spectrum-spreading-signal demodulation apparatus accordingto an embodiment to which a frequency analysis method according to thepresent invention is applied will be described below.

FIG. 14 is a block diagram showing an example structure of the GPSreceiver in that case. As shown in FIG. 14, a signal (spectrum spreadingsignal) sent from a GPS satellite and received by an antenna 21 is sentto an intermediate-frequency conversion circuit 23 through ahigh-frequency amplification circuit 22. The output of a referenceoscillator 24 formed of a crystal oscillator is sent to a localoscillation circuit 25. With this, a local-oscillation output isobtained in which the output frequency and the frequency ratio of thereference oscillator are fixed.

The local-oscillation output is sent to the intermediate-frequencyconversion circuit 23, and a satellite signal is low-frequency-convertedto an intermediate-frequency signal having an intermediate frequency of1.023 MHz. The intermediate-frequency signal is amplified by anamplification circuit 26, its bandwidth is limited by a band-pass filter27, and the signal is sent to a DSP (digital signal processor) 100.

A block diagram enclosed by a dotted line in FIG. 14 shows a functionexecuted by the DSP 100 in a hardware manner. These blocks can also bestructured by discrete circuits as hardware. The structure of the DSP100 in FIG. 14 shows the structure of a digital matched filter.

The signal sent to the DSP 100 is converted to a digital signal by anA/D converter 101, and written into a buffer memory 102. The signal readfrom the buffer memory 102 is sent to an FFT processing section 103. Inthis case, the FFT processing section 103 applies FFT processing to thedigital signal for one period (corresponding to 1023 chips) of aspreading code with the phase being shifted sequentially by one chip.

As the FFT processing section 103, the frequency analysis apparatusaccording to the first embodiment or the second embodiment of thepresent invention, described before, is used. The spectrum spreadingsignal, which is also a received signal in the GPS receiver, is a realsignal, and the second embodiment can be applied as is.

In the GPS receiver in this case, the results of FFT stored in the RAMsof the latter-stage calculation means 2 or 12 are all written into amemory 104. The results of FFT of the received signal, read from thememory 104 are sent to a multiplication section 105.

A spreading-code generation section 106 generates a spreading code ofthe same sequence as the spreading code used in the received signal fromthe satellite, which is to be processed by the DSP 100. The spreadingcode for one period (1023 chips) is sent from the spreading-codegeneration section 106 to an FFT processing section 107.

As the FFT processing section 107, the frequency analysis apparatusaccording to an embodiment of the present invention, described above maybe used. The results of FFT processing are sent to a memory 108. Theresults of FFT are sequentially read from a lower frequency as in ausual case and sent to the multiplication section 105.

The multiplier 105 multiplies the result of FFT processing of thereceived signal sent from the memory 104 by the result of FFT processingof the spreading code sent from the memory 108 to calculate the degreeof correlation between the received signal and the spreading code in afrequency domain. Actually, in the multiplication section 105, thecomplex conjugate of one of the result of discrete Fourier transform ofthe received signal and the result of discrete Fourier transform of thespreading code is multiplied by the other. Specifically, the complexconjugate of the result of FFT read from the memory 108 (or the memory104) is calculated by a calculation section (not shown), and is sent tothe multiplication section 105. Alternatively, before the result of FFTis written into the memory 108 (or the memory 104), the complexconjugate thereof may be calculated. The result of multiplication issent to an inverse-FFT processing section 109, and the signal in thefrequency domain is converted back to a signal in the time domain.

The result of inverse-FFT processing obtained from the inverse-FFTprocessing section 109 is a correlation detection signal of the receivedsignal and the spreading code in the time domain. This correlationdetection signal is sent to a correlation-point detection section 110.The correlation-point detection section 110 determines whether thereceived signal has been synchronized with the spreading code. When itis determined that they have been synchronized, the phase of a peak isdetected as a correlation point.

This correlation detection signal indicates a correlation value at eachchip phase in one period of the spreading code. When the spreading codein the received signal is synchronized with the spreading code sent fromthe spreading-code generation section 106, a correlation waveform isobtained as shown in FIG. 15, in which a correlation value at the phaseof one chip among 1023 chips shows a peak which exceeds a thresholddetermined in advance. The chip phase having the peak is the phase of acorrelation point.

When the spreading code in the received signal is not synchronized withthe spreading code sent from the spreading-code generation section 106,a correlation waveform having a peak as shown in FIG. 15 is notobtained. There is no peak which exceeds the threshold determined inadvance at any chip phase.

The correlation-point detection section 110 determines, for example,whether a correlation detection signal sent to the correlation-pointdetection section 110 has a peak which exceeds the threshold determinedin advance, to determine whether the received signal has beensynchronized with the spreading code. When it is determined that theyhave been synchronized, the phase of the peak is detected as acorrelation point.

In the above description, the carrier of the received signal is notconsidered. Actually, the received signal r(n) includes a carrier asindicated by expression (6) of FIG. 20. In expression (6), “A” indicatesan amplitude, d(n) indicates data, f₀ indicates a carrier angularfrequency in an intermediate-frequency signal, and n(n) indicates noise.

When a sampling frequency is called fs and the number of sampling timesis called “N” (therefore 0≦n<N and 0≦k<N) in the A/D converter 101, therelationship between a discrete frequency k and an actual frequency fafter discrete Fourier transform satisfies: f=k·fs/N at 0≦k≦N/2, andf=(k−N)·fs/N at N/2<k<N. Due to a characteristic of discrete Fouriertransform, R(k) and C(k) are circulative when K<0 and K≦N.

To obtain data d(n) from the received signal r(n), it is necessary tosynchronize a spreading code c(n) with a carrier cos 2πnf₀ and to removethe carrier component. In other words, when only R(k) includes thecarrier component in expression (5) of FIG. 20, described later, awaveform like that shown in FIG. 15 is not obtained.

In the present embodiment, a simple structure just having FFT processingin the frequency domain allows the spreading code c(n) and the carriercos 2πnf₀ to be synchronized and the carrier component to be removed.

Specifically, the results of FFT of a received signal from a GPSsatellite, obtained by the FFT processing section 103, are read from thememory 104 in an ascending order of frequency of the frequencycomponents of the received signal and sent to the multiplication section105 in usual cases. In the present embodiment, however, the results ofFFT of the received signal are sequentially read from the memory withthe reading addresses being shifted by the control of a reading-addresscontrol section 111.

Information on the carrier frequency of the received signal, detectedaccording to a correctly estimated amount of Doppler shift for the GPSsatellite from which the received signal was sent and according to acorrectly calibrated oscillation frequency and time information insidethe GPS receiver is sent to the reading-address control section 111. Theinformation on the carrier frequency can be generated only by the GPSreceiver, but usually is obtained from the outside.

Then, the reading-address control section 111 shifts the readingaddresses by the carrier frequency according to the obtained informationon the carrier frequency, sequentially reads the results of FFT of thereceived signal from the memory 104, and sends them to themultiplication section 105.

When the results of FFT of the received signal r(n) is read from thememory 104 with the reading addresses being shift by the carrierfrequency of the received signal in this way, the results of FFT equalto the results of FFT of the received signal from which the carriercomponent has been removed is obtained, as described later. When theresults of multiplications of the results of FFT of the received signalfrom which the carrier component has been removed and the results of FFTof the spreading code for one period thereof are de-spread, acorrelation detection output having a peak at a correlation point asthat shown in FIG. 15 is always obtained.

As described later, it is also possible that the carrier component ofthe received signal r(n) is added to the results of FFT of the spreadingcode by controlling not the reading address for the results of FFT readfrom the memory 104 but the reading address for the results of FFT ofthe spreading code read from the memory 108, and the carrier componentis substantially removed by the multiplications in the multiplicationsection 105.

The removal of the carrier component enabled by the synchronization ofthe carrier of the received signal and the spreading code by the controlof the reading addresses of the memory 104 or 108 will be describedbelow in more detail together with the operation of the processing of adigital matched filter in the DSP 100.

In the present embodiment, the processing of a digital matched filter isperformed in the DSP 100. The principle of the processing of the digitalmatched filter is based on a theorem in which convolutional Fouriertransform in the time domain is a multiplication in the frequencydomain, as shown in expression (4) of FIG. 20.

In expression (4), r(n) indicates a received signal in the time domain,and R(k) indicates the discrete Fourier transform thereof. In addition,c(n) indicates a spreading code sent from the spreading-code generationsection, C(K) indicates the discrete Fourier transform thereof, “n”indicates a discrete time, “k” indicates a discrete frequency, and F[ ]indicates Fourier transform.

When a correlation function of the two signals r(n) and c(n) is definedas f(n), F(k) which shows the discrete Fourier transform of f(n) has arelationship shown in expression (5) of FIG. 20. Therefore, when it isassumed that r(n) is a signal sent from the A/D converter 101 shown inFIG. 14 and c(n) is a spreading code sent from the spreading-codegeneration section 106, the correlation function f(n) of r(n) and c(n)can be calculated in the following procedure by using expression (5)without using a usual definition expression.

-   -   Calculate R(k) which is the discrete Fourier transform of the        received signal r(n).    -   Calculate the complex conjugate C (k) of C(k) which is the        discrete Fourier transform of the spreading code c(n).    -   Calculate F(k) in expression (2) by using R(k) and the complex        conjugate C (k) of C(k).    -   Calculate a correlation function f(n) by applying inverse        discrete Fourier transform to F(k).

When the spreading code included in the received signal r(n) matches thespreading code c(n) sent from the spreading-code generation section 106,as described above, the correlation function f(n) calculated accordingto the above procedure has a time waveform that has a peak at acorrelation point as shown in FIG. 15.

Next, the synchronization of the carrier included in the received signalr(n) and the spreading code will be described.

As described before, the received signal r(n) includes the carrier asindicated by expression (6) of FIG. 20. To obtain the data d(n) from thereceived signal r(n), it is necessary to synchronize the spreading codec(n) with the carrier cos2πnf₀ and to remove the carrier component. Inother words, when only R(k) includes the carrier component in expression(5) of FIG. 20, described before, a waveform like that shown in FIG. 15is not obtained.

When the amount of Doppler shift is correctly estimated and theoscillation frequency and the time information of the GPS receiver arecorrect, the carrier frequency f₀ of the received signal r(n) becomesknown. In this case, as shown in FIG. 16, it is possible to remove thecarrier component from the received signal r(n) before FFT by providinga multiplication section 121 at a stage prior to the FFT processingsection 103 and by multiplying the received signal r(n) by the carrierhaving the frequency f₀ sent from a signal generation section 122 in themultiplication section 121 for frequency conversion.

In this case, since the results of FFT of the received signal r(n) fromwhich the carrier component has been removed are obtained from thememory 104, and the results of FFT are multiplied by the results of FFTof the spreading code c(n) in the multiplication section 105, a timewaveform having a peak at a correlation point as in FIG. 15 is alwaysobtained as the output of the inverse-FFT processing section 109.

As shown by parentheses in FIG. 16, the same result is obtained when thecarrier component is added to the spreading code without removing thecarrier component from the received signal r(n), by providing themultiplication section 121 at a stage prior to the FFT processingsection 107 and by multiplying the spreading code c(n) by the carrierhaving the frequency f₀ sent from the signal generation section 122 inthe multiplication section 121 for frequency conversion.

In this case, since the carrier component included in the results of FFTof the received signal read from the memory 104 is synchronized with theadded carrier component included in the results of FFT of the spreadingcode read from the memory 108, a correlation detection output having apeak at a correlation point as in FIG. 15 is obtained from theinverse-FFT processing section 109.

In the method described above in which the signal in the time domain ismultiplied by the carrier-frequency signal as shown in FIG. 16, however,the multiplication section used for removing the carrier component isespecially required, the structure becomes complicated, and theprocessing speed becomes slower by the multiplications.

As a characteristic of FFT, the above-described frequencymultiplications can be expressed by expression (7) of FIG. 20. Inexpression (7), F[ ] indicates discrete Fourier transform, φ₀ is a phasedifference from the carrier, k₀ is k corresponding to f₀, andf₀=k₀·fs/N. From expression (7), a signal obtained by applying FFTprocessing to the signal frequency-converted from the received signalr(n) as shown in FIG. 16 is R(k) obtained by applying FFT processing tor(n) and shifted by the carrier frequency k₀.

From the above description, the structure shown in FIG. 16 can bechanged to a structure shown in FIG. 17. More specifically, instead ofmultiplying the received signal r(n) or the spreading code c(n) by thecarrier frequency, the reading addresses used when the results of FFT ofthe received signal or the results of FFT of the spreading code are readfrom the memory 104 or the memory 108 are shifted by the carrierfrequency.

In this case, when the received signal r(n) is shifted, down conversionis used with k₀>0, and when the spreading code c(n) is shifted, upconversion is used with k₀<0, in FIG. 17.

As described above, when the FFT characteristic shown in expression (7)is used, the signal generation section 122 shown in FIG. 16 is notrequired. In addition, as shown in FIG. 17, since only the readingaddress phases used for reading the results of FFT from the memory needto be shifted, the structure is simplified, and the processing is madefaster.

Since the phase difference φ₀ in expression (7) is unknown, it isignored in FIG. 17. For example, a correlation function f′(n) (0≦n<N)obtained as the result of inverse-FFT calculations of F′(k) calculatedby expression (8) of FIG. 20 is a complex number. When its real part iscalled f_(R)′(n) and its imaginary part is called f_(I)′(n), theamplitude |f′(n)| of a correlation peak is obtained as shown inexpression (9) of FIG. 20 and the phase φ is obtained as shown inexpression (10) of FIG. 20. Therefore, the multiplication of exp(jφ₀) atthe right side of expression (7) may be omitted. The phase φ is obtainedby adding φ₀ in expression (7) to two values apart by π, correspondingto the sign of the data d(n) in expression (6)

FIG. 18 shows a structural view obtained when the processing operationdescribed above is applied to the block view of FIG. 14. At the outputsof blocks shown in FIG. 18, signal outputs r(n) and c(n), and theresults R(k), C(k), and f′(n) of calculations like those described aboveare shown.

As described above, according to the demodulation method for spectrumspreading code, having the structure shown in FIG. 14, when FFT is usedto form a digital matched filter in a GPS receiver, the results of FFTof a received signal are multiplied with the spreading code with theaddresses of a memory being shifted by the value corresponding to acarrier frequency, as shown in FIG. 18. A correlation point np isobtained, for example, in a waveform like that shown in FIG. 18. When acorrelation point np is obtained for each of four GPS satellites, thatis, for four spreading codes c(n), the position of the GPS receiver canbe calculated.

In other words, according to the embodiment shown in FIG. 14, when theprocessing of a digital matched filter which uses FFT is performed, thecarrier component of a received signal can be removed by a simple methodin which one of the results of FFT of the received signal and theresults of FFT of the spreading code is shifted when the results of FFTof the received signal and the results of FFT of the spreading code aremultiplied in the frequency domain, without performing multiplicationsin the time domain to synchronize the carrier of the received signalwith the spreading code.

As described before, when the frequency analysis method according to thefirst and second embodiments of the present invention is applied to theFFT processing section 103 and the FFT processing section 107 in thedemodulation apparatus for spectrum spreading signals shown in FIG. 14,the structure of the DSP 100 is simplified, and high-speed processing isexpected. An input signal to the demodulation apparatus for spectrumspreading signals is a real signal, and the second embodiment can alsobe applied as is.

In the above-described cases, the reading addresses of the memory usedfor R(k), which is the results of FFT of the received signal, areshifted. The reading addresses of the memory used for C(k), which is theresults of FFT of the spreading signal, may be shifted (corresponding toup conversion in the multiplication section) in the direction oppositethat used for R(k), which is the results of FFT of the received signal.

In the above description, the spreading-code generation section 106 andthe FFT processing section 107 are separately provided. FFT calculationsfor the spreading code c(n) performed when a satellite signal isreceived can be omitted if the spreading code corresponding to each GPSsatellite is FFT-processed and stored in a memory.

In the above-described embodiment, the carrier frequency of the receivedsignal from the GPS satellite is known. When the carrier frequency isunknown, the same result is obtained by sending the correlationdetection output of the correlation-point detection section 110 to thereading-address control section 11 and by searching for the amount ofshift which is applied to the reading addresses of the results of FFTfrom the memory 104 and which causes a peak like that shown in FIG. 15to be obtained.

More specifically, in this case, the reading-address control section 111changes the amount of shift applied to the reading address of theresults of FFT of the received signal r(n) from the memory 104,according to the correlation detection output of the correlation-pointdetection section 110, around an estimated address determined from pastdata, so that the correlation-point detection section 110 obtains a peaklike that shown in FIG. 15. When the correlation-point detection section110 obtains a peak like that shown in FIG. 15, the reading-addresscontrol section 111 stops the shift control of the reading addresses atthe amount of shift used at that time.

A frequency analysis method according to the present invention can beapplied not only to the above-described GPS receiver but also to variousfrequency analyses.

As described above, according to a frequency analysis method of thepresent invention, calculations for discrete Fourier transform can beperformed at a high speed with a smaller number of memories or a smallercapacity of a memory than in usual FFT, and frequency components can beanalyzed without reducing the frequency resolution.

In this case, a larger amount of calculations is performed than in usualFFT as the memory is reduced, but that amount is greatly smaller thanthe amount of calculations required when the calculations are performedaccording to the definition of discrete Fourier transform.

1. A frequency analysis method using a-power-of-two N_(A)=2^(a), where“a” is an integer, memories for discrete Fourier transform, comprising:a former-stage calculation step of taking out frequency components of aninput signal having a-power-of-two N_(B)=2^(b), where “b” is an integerand b>a, data items, in a comb manner and of calculating N_(A)intermediate data items, a latter-stage calculation step of applyingfast Fourier transform to the intermediate data items obtained in theformer-stage calculation step, by using the N_(A) memories for discreteFourier transform, and executing the former-stage calculation step andthe latter-stage calculation step a predetermined number of times withthe frequency components taken out in the comb manner being changed, toexecute discrete Fourier transform through the predetermined number oftimes for determining the position of a receiver and time information ofa GPS satellite.
 2. The frequency analysis method according to claim 1,wherein the predetermined number of times is 2^(b−a).
 3. The frequencyanalysis method according to claim 1, wherein the input signal is asignal having no imaginary part and in the latter-stage calculationstep, as the results of fast Fourier transform of some data items amongthe intermediate data items, a part of other intermediate data items, inwhich signs of imaginary parts are inverted and data items arere-arranged is output; and wherein processing in the former-stagecalculation step and FFT calculations in the latter-stage calculationstep are omitted for the part of the intermediate data items.
 4. Aspectrum-spreading-signal demodulation method comprising the steps of:applying discrete Fourier transform to a received signal in which acarrier wave is modulated by a signal obtained by spectrum-spreadingdata with a spreading code; detecting correlation between the receivedsignal and the spreading code by multiplying the complex conjugate ofone of the results of discrete Fourier transform of the received signaland the results of discrete Fourier transform of the spreading code bythe other; detecting a correlation point of the received signal and thespreading code by applying inverse Fourier transform to results ofmultiplying, wherein in the step of applying discrete Fourier transformto the received signal, using a-power-of-two N_(A)=2^(a), where “a” isan integer, memories for discrete Fourier transform in executing aformer-stage calculation step of taking out frequency components of aninput signal having a-power-of-two N_(B)=2^(b), where b is an integerand b>a, data items, in a comb manner and of calculating N_(A)intermediate data items, and executing a latter-stage calculation stepof applying fast Fourier transform to the intermediate data itemsobtained in the former-stage calculation step, by using the N_(A)memories for discrete Fourier transform, wherein the former-stagecalculation stip and the latter-stage calculation step are executed apredetermined number of times with the frequency components taken out inthe comb manner being changed, thereby to execute discrete Fouriertransform through the predetermined number of times for demodulating thereceived spectrum-spreading signal.
 5. The spectrum-spreading-signaldemodulation method according to claim 4, wherein the predeterminednumber of times is 2^(b−a).
 6. A spectrum-spreading-signal demodulationmethod according to claim 4, wherein in the latter-stage calculationstep, as the results of fast Fourier transform of a part of theintermediate data items, the results of fast Fourier transform of otherintermediate data items, in which the signs of imaginary parts areinverted and data items are re-arranged, are output; and processing inthe former-stage calculation step and FFT calculations in thelatter-stage calculation step are omitted for the part of theintermediate data items.
 7. A frequency analysis apparatus comprising:a-power-of-two N_(A)=2^(a), where “a” is an integer, memories fordiscrete Fourier transform; former-stage calculation means for takingout frequency components of an input signal having a-power-of-twoN_(B)=2^(b) data items, where “b” is an integer and b>a, in a combmanner and for calculating N_(A) intermediate data items; latter-stagecalculation means for applying fast Fourier transform to theintermediate data items obtained by the former-stage calculation meansby using the N_(A) memories for discrete Fourier transform; andcalculation control means for applying a series of calculationsperformed by the former-stage calculation means and the latter-stagecalculation means to the input signal a predetermined number of timeswith the frequency components taken out in the comb manner beingchanged, thereby to execute discrete Fourier transform through thepredetermined number of times for determining the position of a receiverand time information of a GPS satellite.
 8. The frequency analysisapparatus according to claim 7, wherein the predetermined number oftimes in the calculation control means is 2^(b−a).
 9. The frequencyanalysis apparatus according to claim 7, wherein the input signal is asignal having no imaginary part and the latter-stage calculation meansoutputs as the results of fast Fourier transform of a part of theintermediate data items, the results of fast Fourier transform of otherintermediate data items in which the signs of imaginary parts areinverted and data items are re-arranged; and the calculation controlmeans controls so as not to apply processing in the former-stagecalculation step and FFT calculations in the latter-stage calculationstep to the part of the intermediate data items.
 10. Aspectrum-spreading-signal demodulation apparatus comprising:received-signal Fourier transform means for applying discrete Fouriertransform to a received signal in which a carrier wave is modulated by asignal obtained by spectrum-spreading data with a spreading code;multiplication means for multiplying the complex conjugate of one of theresults of discrete Fourier transform of the received signal and theresults of discrete Fourier transform of the spreading code used in thereceived signal, by the other; inverse Fourier transform means forapplying inverse Fourier transform to the results of multiplicationsobtained by the multiplication means to obtain a correlation detectionoutput of the received signal and the spreading code; and means fordetecting a correlation point of the received signal and the spreadingcode by searching for a peak of correlation of the received signal andthe spreading code according to the correlation detection outputobtained by the inverse Fourier transform means, wherein thereceived-signal Fourier transform means includes: a-power-of-twoN_(A)=2^(a), where “a” is an integer, memories for discrete Fouriertransform; former-stage calculation means for taking out frequencycomponents of an input signal having a-power-of-two N_(B)=2^(b) dataitems, where b is an integer and b>a, in a comb manner and forcalculating N_(A) intermediate data items; latter-stage calculationmeans for applying fast Fourier transform to the intermediate data itemsobtained by the former-stage calculation means, by using the N_(A)memories for discrete Fourier transform; and calculation control meansfor applying a series of calculations performed by the former-stagecalculation means and the latter-stage calculation means to the inputsignal a predetermined number of times with the frequency componentstaken out in the comb manner being changed, to execute discrete Fouriertransform through the predetermined number of operations fordemodulating the received spectrum-spreading signal.
 11. Thespectrum-spreading-signal demodulation apparatus according to claim 10,wherein the predetermined number of times in the calculation controlmeans is 2^(b−a).
 12. The spectrum-spreading-signal demodulationapparatus according to claim 10, wherein the latter-stage calculationmeans outputs, as the results of fast Fourier transform of a part of theintermediate data items, the results of fast Fourier transform of otherintermediate data items, in which the signs of imaginary parts areinverted and data items are re-arranged; and the calculation controlmeans controls so as not to apply processing in the former-stagecalculation step and FFT calculations in the latter-stage calculationstep to the part of the intermediate data items.
 13. Thespectrum-spreading-signal demodulation apparatus according to claim 10,wherein the received-signal Fourier transform means is formed of adigital signal processor.